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Abstract 

We compare some first order well-balanced numerical schemes for shallow water system with 
special interest in applications where there are abrupt variations of the topography. We show 
that the space step required to obtain a prescribed error depends on the method. Moreover, 
the solutions given by the numerical scheme can be significantly difi'erent if not enough space 
resolution is used. We shall pay special attention to the well-known hydrostatic reconstruction 
technique where it is shown that large bottom discontinuities may be neglected and a modification 
is proposed to avoid this problem. 



1 Introduction 

The numerical analysis of the shallow water system has been extensively studied in recent years. This 
system may be written in the form: 

dh dq^ 
dt dx ' 

(1.1) 

The variable x makes reference to the axis of the channel and t is time; q{x, t) and t) represent 
the mass-flow and the thickness, respectively; 5, the acceleration due to gravity; H{x), the depth 
measured from a fixed level of reference; q{x,t) = h{x,t)u{x,t), with u the depth averaged horizontal 
velocity. 

This model has the structure of a system of balance laws: 

dtw + d^F{w) = Siw)d^H, (1.2) 

with 

By considering H as an (artificial) unknown whose value is given by the initial conditions, the system 
can also be written as a first order quasilinear system: 

Wt+A{W)-W^^ = 0, xeR,t>0, (1.3) 
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with W = {h, q, H)* and 



where 



/ 1 \ 

A{W) = -u^ + 2m -c2 , (1.4) 
\ / 



w=|, c=y/gh. 
The eigenvalues of ^(W) are 

Ai = M — c, A2 = M + c, A3 = 0, 

The system is strictly hyperbolic provided that h > Q and |m| 7^ c. 
The Riemann invariants corresponding to the null eigenvalues are: 

q = constant, h + „ — H = constant. (1.5) 
2gh'' 

In the particular case q = these curves are straight lines in h, q, H space: 

g = 0, h — H = constant, (1.6) 

Once the depth function H{x) has been fixed, the stationary solutions of the system are given by 
(1.5). In particular, those given by (1.6) correspond to water at rest solutions. 

Models of the form (1.2) or (1.3) also appear in many other physical problems. It is well known that 
standard methods that solve correctly systems of conservation laws can fail in the presence of source 
terms, especially when approaching equilibria or near to equilibria solutions. In the context of shallow 
water equations Bermiidez and Vazquez-Cendon introduced in [24] the condition called conservation 
property or C-property; a scheme is said to satisfy this condition if it solves correctly the steady state 
solutions corresponding to water at rest. This idea of constructing numerical schemes that preserve 
some equilibria, which are called in general well-balanced schemes, has been studied by many authors 
(see, for instance, [14], [6], [19] [11], [17], [2], [7], [3], among many others). 

This study was motivated by the following observation: given a fixed mesh, diS'erent well-balanced 
numerical methods based on different treatments of the source terms can lead to very different nu- 
merical solutions. Figure 1 illustrates this phenomenon: a numerical test consisting of a shallow layer 
of fluid flowing up a ramp has been considered, and the water elevation obtained with different first 
order numerical schemes using the same mesh are compared. Some similar test cases as well as the 
considered numerical schemes will be described in detail below. As it can be seen, very different be- 
haviors of the fluid layer is predicted by the considered numerical schemes; in some cases the fluid goes 
up the ramp in a smooth way, in some others a strong shock going to the left is predicted. It is worth 
noting that these huge differences are not related to a convergence problem: these differences can be 
done arbitrarily small by refining the mesh and/or by increasing the order of accuracy of the method. 
Nevertheless, for different first order methods, the refinement level required to obtain a solution close 
to the common limit can be significantly bigger than those necessary for some others. Moreover, when 
solving these problems with a first order finite volume method, one approaches the solution with a 
piece-wise constant function. This means that in practice we have a discontinuity of the bottom at 
each interface, which motivates the study not only for steep topography but in the presence of a 
discontinuous bottom, where we shall see that this kind of phenomenon arises. This phenomenon 
implies a difficulty: in practical applications one usually has a given mesh whose refinement level 
depends on the available bathymetric data and on the available computational power. Therefore, it 
is important to know how reliable the numerical solutions are if there are abrupt variations of the 
topography and the mesh is not fine enough. Roughly speaking, a numerical treatment of the source 
term is expected to be more reliable than others if the refinement level of the mesh which is required 
to obtain a numerical solution whose behavior is close enough to the limit is lower. 
The goal of this article is to run a battery of numerical tests in order to compare different numerical 
schemes based on different treatment of the source terms and to draw conclusions about their relia- 
bility. The different numerical treatment of the source terms considered here belong to two families: 
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(a) Surface profiles (b) Surface profiles 

Figure 1: Comparison of numerical solutions obtained with Roe scheme and hydrostatic reconstruc- 
tion. Between parenthesis the number of cells used is shown 



on the one hand, we consider numerical methods based on an upwinding treatment of the source 
term. This treatment has been introduced by Roe for scalar problems [20] and by Bermiidez and 
Vazquez-Cendon in [24]. On the other hand, we consider numerical methods based on the hydrostatic 
reconstruction technique introduced in [2] and some extensions and modifications. A phenomenon 
which is closely related to the one that motivated this paper was observed in [13] for numerical 
methods based on the hydrostatic reconstruction technique: it was seen that the numerical solutions 
obtained on a fixed mesh with a first order numerical scheme for a test consisting of a shallow layer 
of fluid going down a ramp were independent of the slope beyond a critical threshold, what is not in 
good agreement with the physics of the problem. This difficulty was also discussed in [5]. A modi- 
fication of the hydrostatic reconstruction technique will be introduced here that partially overcomes 
this difficulty. 

Although the immerical analysis of system (1.1) with discontinuous depth function H is not addressed 
here, the theoretical and numerical difficulties arising in this case can give some insight about the 
origin of the phenomenon to be studied here. When the depth function H is discontinuous, the source 
term cannot be defined within the distributional framework: it is a non-conservative product that 
may be defined in many different ways: see, for instance, [12]. Depending on the chosen definition 
of the source term across a discontinuity of H, one has many different ways to define what is a weak 
solution of the problem. According to the more natural and extended definition, the simple waves 
associated to a bottom discontinuity are contact discontinuity preserving the Riemann invariants of 
the third characteristic field (1.5): see [14], [1], [15]... Nevertheless, some alternative definitions of 
weak solution can be found in the literature: see [4], [21]... 

Once the notion of weak solution has been fitted, a concept of entropy solution is necessary to discard 
the unphysical ones. In the case of the shallow water system, there is a known entropy pair that allows 
one to define such a concept: a weak solution is said to be an entropy solution if the inequality 

dtv{w,H) + d,G{w,H)<0 (1.7) 

is satisfied in the distribution sense, where 

VH = h— + G{w) =[j+ ghj hu, 

ri{w, H) = r]{w) - ghH, G{w, H) = G{w) - ghuH. 

Nevertheless, in this case the addition of an entropy notion does not ensure the uniqueness of solution: 
in the case of resonance problems, i.e. when one of the eigenvalues of the matrix system changes its 
sign: see, for instance, [15]. 
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From the numerical point of view, the main difficulty of solving the system with a discontinuous 
depth function comes from the fact that the solutions provided by two different numerical methods 

can converge to functions which are weak solutions according to different definitions, or to different 
entropy solutions according to the same definition of weak solutions. 

All the first order numerical methods compared in this article are based on piecewise constant dis- 
cretizations of the depth function H. Therefore, the solution of a Riemann problem related to the 
non-conservative system (1.3) with discontinuous values of H has to be approached at every intercell. 

As different approximate Riemann solvers may converge to different solutions, these local approxi- 
mations may be very difi^erent to each other: this is the reason of the different behaviors observed in 
Figure 1. 

This paper is organized as follows: in Section 2 and for the sake of completeness, the different numerical 
schemes to be compared are briefly presented. Section 3 is devoted to the hydrostatic reconstruction, 
focusing on the difficulty mentioned before and a modified version of the hydrostatic reconstruction 
technique is proposed to overcome it. Finally, Section 4 will present different numerical simulations 
comparing previously cited numerical schemes. 



2 Upwind treatment of the source term 



We consider as usual a set of computing cells Ii = [.Xj_i/25 a;i+i/2]) * G ^- We shall assume, for the 
sake of simplicity, that these cells have a constant size Ax and that = iAx. Xi = {i — l/2)Ax 

is the center of the cell 7j. Let At be the time step and = nAt. 

The following notation will be used for the approximation of the cell averages of the exact solution: 



wp 



Hi 



We consider numerical schemes that can be written in the form 



At 



At 



where 



(2.1) 

(2.2) 
(2.3) 

(2.4) 
(2.5) 



F,+l/2=^«,<+l), 

being a continuous function such that 

J^{w,w) = F{w), 

and 

where S~ and 5+ have two components: a centered approximation of the source term and the upwind 

treatment. 

Let us introduce some examples: 
2.1 Roe method 

The extension of the Roe method to (1.2) is given by (2.2) with the choices 



where 



'^i-l-1/2 - 'Pi+l/2^i+l/2iHi+l - Hi), 

1 



Ji+1/2 



(2.6) 
(2.7) 

(2.8) 
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where 



Si+\/2 is given by: 



and 




Here Id represents the identity matrix and 

\Ji+ll2\ = ■'fi+l/2|Aj+l/2|-ffj^\/2' (2-12) 

where |Aj_|_i/2l is the diagonal matrix wliose coefficients are the absolute value of the eigenvalues of 
Jj+i/2 and Ki^ii2 is the matrix whose columns are associated eigenvectors. This numerical scheme 
is exactly well-balanced for water-at-rest solutions and is second order accurate for general stationary 
solutions: see [24], [19]. 

2.2 FORCE and GFORCE methods 

FORCE and GFORCE schemes were introduced in [22] and [23] for conservative systems. Their well- 
balanced extensions to nonconservative systems were introduced in [10] . For the particular case of the 
shallow water system, these extensions can be written in the form (1.2) with 



j;+l/2=;.(i^«)+J^«+l)) 



Ax . 
At 



At . 
Id + i^—Jt 



(2.13) 



Si+l/2{Hi+l — Hi, 



, , ,Ax At ^ 

' (^"'^^ At *+i/2+'^ Ax 



Si+,/2{Hi+i - Hi) (2.14) 



where Ji+1/2 and S'j_|_i/2 are given by (2.8) and (2.10) respectively. FORCE and GFORCE methods 
correspond, respectively, to w = 1/2 and w = 1/(1 + CFL), where CFL G (0, 1] is the CFL number. 
Notice that the choices ui = and w = 1 correspond to extensions of the Lax-Friedrichs and Lax- 
Wendroff methods. This numerical scheme is also exactly well-balanced for water at rest solution and 
second order accurate for general stationary solutions. 

Remark that (2.14) is not well defined when one of the eigenvalues of Ji+1/2 vanishes. This means 
that when there is a change from subcritical to supercritical regime, this scheme will give problems. 
To correct this, we may replace 



„2\ I7T„2 _ „,2 



)2 
2> 



2u. 



'i-H/2 



(2.15) 



where /i = max{e, jl — Fr^j} ■ sgn(l — Fr^), Fr^ = M?|_j/2/(ci-i-i/2)^ is the square of Froude number 
and e > is some given small value. 

Another possibility would be to use the following modification of the discrete source term: 



'^i-l-1/2 



'S'i+l/2(-ffi+l — Hi 



±((l-a;)^(J^i/2) 



At ^ 



Si+^/2{Hi+r~Hi) (2.16) 
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where: 



7* 



1/2 







This modified method is still exactly well-balanced for water-at-rest solutions but only first order 
accurate for general stationary solutions. We refer to [10] for further details. 



2.3 Path-conservativity 

The upwind numerical methods introduced in this section satisfy the following equalities: 

Sf^^l^=QifWr = WJX„ (2.18) 

and 

5+1/2 + 5-^/2 = f^S{^^{s-Wr,Wr+r))^[s-Wr,W:;Xi)ds, (2.19) 

where 

^(.s: H'7, W'^_^-^) = WP + s {WlXi - WP) , (2.20) 

with ^I^TO and '^h the two first and last component of * respectively. These properties imply that 
these methods can be interpreted as path-conservative numerical methods in the sense defined in [18] 
for (1.3), that is, these methods are formally consistent with a particular definition of weak solution of 
the Riemann problem with a discontinuous depth function H: more precisely, they are consistent with 
the definition of nonconservative product provided by the theory developed in [12] corresponding to 
the choice of the family of straight segments (2.20). The choice of the family of straight segments is a 
natural choice to obtain well-balanced numerical methods for water-at-rest solutions, as the equation 
(1.6) corresponds to straight lines in the h,q,H space: [16] for details. 



3 Hydrostatic reconstruction technique 
3.1 The original formulation 

The hydrostatic reconstruction technique was introduced in [2] . It presents a straightforward method 
for obtaining a well-balanced scheme for the shallow water system based on a numerical flux for the 
homogeneous system H = constant. Moreover, this technique has the advantage of preserving good 
properties of the homogeneous solver such as positivity of water thickness and entropy inequalities. 
The idea is to consider a consistent numerical flux J^{wi,Wr) for the homogeneous shallow water 
system, that is (1.1) with H constant, and the corresponding numerical scheme 

<+l = < + ^{Fi-l/2 - Ji+l/2), Fi+y2 = Hwi, Wi+l). (3.1) 

Then define at each intercell the values 

-ffi+i/2 = min(ffi,i?i+i/2), 1i+i/2 = K+\/2'^i^ 9i+i/2 = '*i+i/2'^'+i' (^-2) 

'*i+l/2 = (hi+i - Hi+\ + _ffj+i/2) + , 
'^i+1/2 = - Hi + -ffi+l/2)+, 

and define the scheme for non-homogeneous shallow water equations in the form (2.2) with 



(3.3) 



Fi+i/2 = •^(Wi+i/2. «',^i/2) (3-4) 

^^^+1/2 = ( P{hi) - p(V^i/2) ) ' ^iv2 ( p{hi+,) -P{h+^y,) ) ' (^-^^ 

where p(h) = ^fi^. 

The resulting numerical method is non-negative and semi-discrete entropy satisfying provided that 
the homogeneous flux also satisfies these properties. 
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3.2 Path-conservativity 



The resulting scheme can be also interpreted as a path-conservative method for (1.3), i.e. (2.18)-(2.19) 
are again satisfied but for a family of paths different from the straight segments (2.20). The family of 
paths in now defined as follows: given WJ^ and W-^^, the path linking them is composed by: 

1. the segment of the straight line of the space h, q, H defined by the equations: 

q = hu'l, h-H = h1-H,, 



linking the states Wp and 



i+l/2 



"i+1/2 
<V+l/2 



-1/2 



2. the segment linking the states W^^^j^ ^^id 



i+l/2 



h+ 

' 'i+l/2 
"i+l"'i+l/2 



H.. 



i+l/2 



3. the segment of the straight line of the space h, q, H defined by the equations: 



"■■i+i 



linking the states W^-^^^ ^^'^ ^P+i- 



Notice that when H, 



i+l/2 



Hi the first segment degenerates to a point and when H, 



i+l/2 



Hi^i is the 



third one that degenerates to a point. Figure 3(a) shows the third segment in the case -ffi+1/2 = Hi. 
The discretized source terms (3.5) are the integral of the source term 



ghdH, (3.6) 



computed by using the first and third segments as the integration path C. Observe that the second 
one does not contribute to the source term as it lies in a plane H = constant. 

When i?i+i/2 = Hi and /i^+i - -ffj+i + Hi < or 7?j+i/2 = -ffj+i and hf - Hi + Hi+i < 0, the paths 
are modified to grant the positivity of the numerical scheme. Let us consider the first case, the second 
one being analogous: the situation is that shown in Figure 2. In this case, h^_^-^/2 ~ ^'^'^ path 

























i 



Figure 2: Case /i+ 



+1/2 



linking the states is composed by: 
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h = {h,+i - H,+i + H)+ 




1 + 1/2 




(b) Case //^_]_i/2 — ^J+i ^ < 0: original 

(dashed line) and modified reconstruction technique (solid 
line) . 



Figure 3: Hydrostatic reconstruction technique. 



1. The segment linking and 



2. the arc of the graph 



hnking the states W 



i+l/2 



h = {hi+i 
and WJ\, . 



i+l/2 



H,,- 






Hi 



Figure 3(b) shows the second piece of this path (dashed hue). 

Observe that, in this case, the only piece of the path that contributes to the discretized source term 
is the segment linking the state W-^i and [0, 0, -ffj+i — ^F+i]"^ in which the straight line 



hu 



i+i 



intersects the axis 



0, 



hu. 



i+l- 



Therefore, the numerical solution does not depend on the value of Hi in this case: this is the reason 
of the phenomenon observed in [13] and mentioned in the Introduction. 

The interpretation of the hydrostatic reconstruction technique in terms of the choice of a particular 
family of paths has been on the basis of the generalization of the technique introduced in [11] that 
allows to obtain a first order numerical method which is exactly well-balanced for every stationary 
solution. It is also on the basis of the simpler modification proposed in the next subsection to partially 
overcome this loss of accuracy in the presence of big steps. 



3.3 Modified hydrostatic reconstruction 

In order to take into account the amplitude of the step, we propose a modification of the hydrostatic 
reconstruction technique. In terms of choice of the family of paths, the proposal is applied whenever 

Hi^i and h" — Hi + Hi^i < 0. In the first case, 



H, 



i+l/2 



Hi and — Hi. 



Hi<0- or H, 



i+l/2 



the idea is to replace the path described above by: 



1. The segment linking and 



i+l/2 






Hi 
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Hi 




Figure 4: Emerging bottom 



2. the segment linking W^j^-^j^ and WJ^j. 

(See the sohd line in Figure 3(b)). The analogous modification is proposed when = fi^i+i and 

/iF - H.. + H,+i < 0. 

In order to keep a similar notation to the one used in the original hydrostatic reconstruction, the 
numerical scheme is written in the form (2.2) where now: 

Fi+1/2 = HWi+l/2:^t+l/2) (S-'') 

s- -( ° 

^++1/2 = [ p(/i,+i) -P(/i++i/2) +7;Xi/2 

Some easy computations of the integral (3.6) through the modified paths lead to the following results: 

1. General case. 

If hi — Hi + min(_ffj, _ffj+i) > and /ij+i — iJj+i + min(_ffj, _ffj+i) > 0, then the standard 
hydrostatic reconstruction is used. Therefore: 

2. Large steps at the intercell. 

If hi — Hi + mm{Hi, Hij^i) < or — -ffj+i + min(_ffj, -ffj+i) < 0, then the path is modified 

= PiK+1/2) - P(h.) + /' ^ (gi+1/2 - H,), (3.9) 

Ti'+i/2=Pihti/2)-Pih'^+i)+g^^^^^^P^iH.+i - H.+-i/2). (3.10) 

The resulting numerical method takes into account the total height of the step and it is well-balanced 
for water at rest solutions, provided that no wet/dry fronts are present, as in this case, the method 
coincides with the original one. Nevertheless, this modification can produce unphysical solutions in 
the presence of an emerging bottom: see Figure 4. 
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Consider for example the case of emerging bottom to the right (Figure 4(a)) and assume that water 
is at rest m = 0. It can be easily checked that this steady state is not preserved by the modified 
hydrostatic reconstruction. On the other hand, if m, > in this situation, two different cases will be 
distinguished: 

• The mechanical energy at the wet cell is enough to make the fluid to go up the step. Then, the 
full step should be considered when computing the numerical fluxes. 

• There is not enough energy at the wet cell, and thus the step acts as an obstacle. 

Following the ideas presented in [8], the fluid should go up the step in the case of emerging bottom 
to the right and Wj > if 

+ g{hi - H, + Hi+i) > IViahiUiY, (3.11) 
and in the case of emerging bottom to the left and Uj+i < if 

+ g{hi+, - Hi+i + Hi) > ^^{ghi+.Ui+^r (3-12) 
Therefore, in the presence of emerging bottom situations, the numerical method is applied as follows: 

• In water at rest situations or if the mechanical energy is not enough to make the water go up 
the step, the original hydrostatic reconstruction technique is applied. 

• If the mechanical energy is enough, the modified reconstruction technique is applied to take into 

account the total height of the step. 

The non-negativity-preserving property of the numerical methods based on the hydrostatic recon- 
struction technique is also satisfied for the modified technique: 

Proposition 3.1. Consider a numerical scheme in the form (2.2). (3. 7). (3. 8) with the modified 
hydrostatic reconstruction described above. If T preserves the non-negativity of h by interface, then 
so does the modified hydrostatic reconstruction. 

The proof of this result is essentially the same as the one presented for hydrostatic reconstruction in 
[6]. 

The entropy preserving property is discussed in Appendix A. 
3.4 Subsonic reconstruction scheme 

The subsonic reconstruction scheme was introduced in [7] as a generalization of the hydrostatic recon- 
struction. This solver intends to preserve a further class of steady state solutions for the shallow water 
system while preserving the good properties of being positive and entropy satisfying. The definition 
of this scheme is rather technical and for the sake of extension it is not included here. We refer the 
reader to [7] for a complete description. 

4 Numerical simulations 

A battery of numerical tests have been designed to compare the numerical solutions obtained with 
the different discretizations of the source terms introduced in Sections 2 and 3. The original or the 
modified reconstruction technique are used by default combined with Roe's flux for the homogeneous 
problem, but in some cases, the FORCE or the GFORCE numerical fluxes are considered. The 
following labels are used in the Figures: 

• Roe: the numerical method given by (2.2), (2.6), (2.7). 

• HR: the hydrostatic reconstruction technique combined with Roe's numerical flux for the homo- 
geneous problem. 
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(a) HR scheme 



(b) Modified HR 



Figure 5: Test 1: water height for 50 points for a = 16, 17, ... , 21 

• Modified HR: the modified hydrostatic reconstruction technique combined with Roe's numerical 
flux for the homogeneous problem. 

• FORCE HR: the hydrostatic reconstruction technique combined with FORCE numerical flux 
for the homogeneous problem. 

• FORCE WB: the numerical method given by (2.2), (2.13), (2.14), and uj = 1/2. 

• GFORCE HR: the hydrostatic reconstruction technique combined with GFORCE munerical flux 
for the homogeneous problem. 

• GFORCE WB: the numerical method given by (2.2), (2.13), (2.14), and cj = 1/(1 + CFL). 

• Subsonic rec: the numerical method described in [7]. 

4.1 Test 1 

This first test intends to show the problem observed in [13]. To do so, let us reproduce the test shown 
there by considering in the interval [0, 3] the initial condition 

h{x,t = 0) =0.02, (j(x,t = 0) = 0.01, H(x) = ^x, (4.1) 

for different values of slope a. Boundary conditions to the left are given by q{x = 0,t) = 0.01, h{x = 
0, t) = 0.02, and free boundary conditions are set on the right boundary. 

As it was remarked by Delestre, water height differs up to a given critical value of the slope. For 
instance, in Figure 5(a) we see that when we use 50 points in the interval [0, 3], HR scheme reproduces 
always the same water thickness. This is related to what has already been discussed before: for the 
resolution of 50 points, we found at each intercell at discrete level a jump which is large enough for 
HR scheme to neglect it. The modified HR scheme corrects this fault as it is shown in Figure 5(b). 
This can be corrected also by using a second order reconstruction technique as it was shown in [13]. 
We remark that here the bottom is continuous, so one expects that when Ax — this problem will 
not be found as we see in Figure 6. 

4.2 Test 2 

This test shows that all the numerical methods give similar result when the depth function is con- 
tinuous and the steps of the piecewise constant discretization of the bottom are small enough. We 
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(a) HR scheme (b) Modified HR 

Figure 6: Test 1: water height for 150 points for a = 16, 17, ... , 21 






— Roe 

f- -T Subsonic rec. 
Modified HR 

H — H HR 

— Bottom 







(a) Surface elevation at time t — 200 



(b) Surface elevation: zoom 



Figure 7: Test 2: transcritical solution 



consider a test taken from [6] consisting of a transcritical flow with a shock over a bump. The interval 
is [0, 25], the initial condition and depth bottom is given by 



Hix) 



h{x, t = 0) = 0.33, q{x, t = 0)= 0.18, 
-0.2 + 0.05(x- 10)2, if8<.i:<12. 



(4.2) 



0, 



otherwise. 



and the boundary conditions are given by q{x = 0,t) = 0.18, h{x = 25, t) = 0.33. 
We consider a uniform mesh composed by 200 cells. The height of the steps of the discretized bottom 
are lower than 1/200. The numerical results are shown in Figure 7. As expected, all the numerical 
results are similar. 



4.3 Test 3 

As it was mentioned in the Introduction, in the presence of a discontinuity of the bottom, the nu- 
merical solutions provided by different numerical methods may converge to different limits which are 
weak solutions of the problem according to different notions, or to different entropy solutions. This 
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— Roe 

T- -T Subsonic rec. 
^ Modified HR 



— Roe 

»- -T Subsonic rec. 
Modified HR 



-T-V -4-1^ -T-T- V-T-'-T-V - 



(a) = 0.15, t = 0.50 (b) Hr = 0.45, t = 3 

Figure 8: Test 3: solutions obtained for two different values of 7?^ in (4.3) 



phenomenon is studied in this Test: we consider the following initial condition and depth function 

h{x,t = 0) = 0.1, q{x,t = 0) =0.15, H{x) = l (4.3) 
^ ' ^ ^ ' i/,,, otherwise, ^ ' 

with increasing values of Hy. from 0.1 to 0.5. The left boundary conditions are given by h{x = 0,t) = 
0.1, q{x = 0,t) = 0.1. and open boundary conditions are considered at x = 1. A uniform mesh of 
200 cells is considered. Figure 8 shows the water elevation obtained for = 0.15 and H,. = 0.45 
at different times: the differences of the numerical solutions provided by the considered munerical 
methods are apparent. 

In order to systematically study these differences, we compute the values of the water thickness to 
the right of the discontinuity obtained with the different numerical methods when a steady state is 
reached. They are also compared with the value h.^ such that the Riemann invariants corresponding to 
the states Wi = (0.1, 0.1, 0.1) and Wr = {hr, 0.1, Hr) coincide: as it was mentioned in the Introduction, 
according to the most used definition of weak solution, the exact stationary solution of the problem 
will be given by a contact discontinuity linking the states Wi and Wr- Therefore, this value of hr is 
labelled as 'exact' in the Figures. 

The results are shown in Figure 9. As it has been reported previously, the results given by the 
hydrostatic reconstruction technique are the same once the height of the step are greater than a 
critical value. The modification introduced here corrects this problem. Nevertheless, we remark that 
none of the schemes reproduce the exact solution, specially in presence of big discontinuities of the 
bottom. This is a known fact that has been discussed in [9] . 

4.4 Test 4 

We consider again a step-like bottom, but in this case the water layer flows from the bottom to the 
top of the step. The initial condition is given again by (4.3) but the depth function is now: 

1 Hr, otherwise , 

with Hr varying from 0.8 to 0.4. A uniform mesh of 200 cells is again considered. The boundary 
conditions h{x = 0,t) = 0.1, q(x = Q,t) = 0.15 are imposed at a; = and open boundary conditions 
are considered at a; = 1. Different munerical solutions are shown in Figures 10 to 12: as it can be 
seen, the bigger the step, the bigger the differences. In order to compare the different behaviors, 
Fr = v? I (gh) is also depicted. 
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Exact 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



(b) 



Figure 9: Test 3: water thickness to the right of the step when a steady-state is reached 
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Figure 10: Test 4: Hr = 0.4, t = 4.5 
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(a) Surface elevation 

Figure 11: Test 4: Hr = 0.65, t = 2 



(b) Fr2 
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Figure 12: Test 4: Hr = 0.75, t = 5) 
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0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 



(a) Left state hi 




0.45 0.5 0.55 0.6 0.65 0.7 0.75 



(b) Right state hr 



Figure 13: Test 4: Comparison of the thickness to the left and to the right of the step once a steady 
state has been reached 



As in the previous test, we compute the values of the thickness hi and to the left and to the right 
of the step once a steady state has been reached. The results are shown in Figure 13 and Figure 14. 
Both hydrostatic reconstructions give essentially the same results but they are very different from the 
ones obtained with the subsonic reconstruction scheme or Roe. Remark that the value hr is the same 
for all the schemes for Hj. = 0.8, which is the case corresponding to a flat bottom. The same occurs 
for Hr < 0.5 when all the numerical solutions predict a reflected shock. 

In this case, we have transcritical regions and (2.15) should be use. We remark that results shown in 
Figure 14(b) depend on the parameter e. Moreover, in this particular case, if we use (2.16), the results 
provided by FORCE WB and GFORCE WB are nearer to those given by their HR version. 

4.5 Test 5 

We consider now a test consisting of a shallow layer of water flowing up a ramp with increasing slope. 
More explicitly, we consider the following depth function and initial condition: 
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0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 



(a) Left state hi 



(b) Right state hr 



Figure 14: Test 4: Comparison of the thickness to the left and to the right of the step once a steady 
state has been reached 



1, if X < xi, 
2 8 

H{x) = { 1 + ^— — (a; - Xi), if < s < 4, (4.5) 
0.2, if a; > 4 

/.(.,. = 0) = (i?-0.9),, ^i^^t = ,) = [l'2J<^^ = '^>'' , (4.6) 

with boundary conditions h{x = Q,t) = 0.1, q(x = 0,t) = 0.9, and xi < A being a given parameter. 
Figure 15 shows the bottom and the initial surface elevation for xi = 3.75. 



- Surface (initial) 

- Bottom 



Figure 15: Test 5: example of initial condition 



Figures 16 to 21 show the computed water thickness before the slope hi and after the slope hp. at 
time t = 2.5 for three uniform meshes of the interval [0, 5] composed by 200, 400, and 800 cells. 
Remark that, as the slope becomes steeper, the values obtained with both the hydrostatic reconstruc- 
tions are further from those obtained for the other schemes. The bigger number of cells, the lesser 
the differences, as one would expect. Moreover, if the grid is not fine enough, extremely different 
behaviors can be observed to the left of the slope; either the water goes up the ramp in a smooth way 
or a shock traveling to the left develops: see Figure 22. 
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(a) hi (b) hr 

Figure 21: Test 5: Number of points: 800 
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(a) Surface elevation (b) Fr^ 

Figure 22: Test 5: xi close to 4, t = 2.5 



4.6 Test 6 

We consider finally a test similar to the previous one, consisting of a shallow layer of water flowing 
down a ramp with increasing slope. More explicitly, we consider the following depth function and 
initial condition 



{0.1, if X < 0.2, 
0.1 + - 0-2), if 0.2 < X < Xr, (4.7) 

Xj- U i ^ 
Hr if X > Xr 

with boundary conditions h{x = 0,t) = 0.5, q{x = 0,t) = 1.2, Xr = 0.2 + Al, Hr = 0.1 + AH with Al 
and AH given parameters. 

The smooth stationary solution reached can be computed from (1.5) and it is shown, together with 
the initial condition, in Figure 23 for a particular value of Al and AH. 

The numerical methods are applied to this test with different values of Al and AH and different 
number of cells in the interval [0, 5]. More explicitly, 8 uniform meshes having 100, 200, 400, 800, 1600, 
3200, 6400, and 12800 cells have been considered. Then, the exact stationary solution is computed for 
every value of Al and AH and it is compared with the ones computed with the numerical schemes. 
As expected, all the numerical stationary solutions converge to the exact one as the mesh is refined. 
Nevertheless, the number of cells required to satisfy a given error bound strongly depends on the 
numerical method: Figures 24 to 26 show the number of cells needed in order to obtain a L^-error 
lower than 0.008. Incomplete graphs mean that the prescribed error bound is not reached for the finer 
grid of 12800 cells. 

As it can be seen, Roe or the subsonic reconstruction schemes need only a small number of cells 
to obtain the prescribed bound, even for the steepest ramp. But this is not the case for the other 
methods: this number dramatically increases as the slope of the ramp increases. 



5 Conclusions 

We have compared different well-balanced numerical schemes for the shallow water system based on 
different discretizations of the source term in test cases in which either the bottom is discontinuous 
or has steep slopes. While the numerical solutions can converge to different limits for discontinuous 
bottom, this is not the case for continuous depth functions. Nevertheless, the space step required 
to obtain an error which is lower than a prescribed bound dramatically depends on the method. 
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Figure 23: Test 6: Surface elevations corresponding to the initial condition and to the stationary 
solution. 



Moreover, if the mesh is not fine enough, the numerical solutions can exhibit completely different 
behaviors. A numerical method is considered to be more reliable as the space step required to achieve 
a prescribed accuracy is bigger. 

We have considered two different approaches to discretize the source term: an upwind treatment or 
the hydrostatic reconstruction technique. While this latter technique, introduced in [2], allows one 
to easily define a robust scheme for shallow water equations, its reliability (in the sense mentioned 
above) is not very high: this is due to the fact that it neglects large bottom discontinuities. We 
have introduced a modification of this technique that takes into account the total height of the 
bottom discontinuities. . Nevertheless, the numerical tests show that both the original and the 
modified hydrostatic reconstruction are less reliable than other numerical methods such as the subsonic 
reconstruction scheme [7] or the Roe scheme [19]. In general, it has been observed that, for continuous 
steep bottoms, a large number of discretization points in the domain is needed for numerical methods 
based on the hydrostatic reconstruction techniques, FORCE or GFORCE methods, while Roe or the 
subsonic reconstruction schemes give good accuracy even with a small number of cells. 

A Appendix 

In [2] it was shown that whenever T is an entropy satisfying numerical flux for the homogeneous shallow 
water system, the hydrostatic reconstruction scheme satisfies a semi-discrete entropy inequality for 
the full system in the sense given by the following definition: 

Definition A.l. A numerical scheme for (1.2) in the form (2.2) - (2.3) is said to satisfy a semi- 
discrete entropy inequality with respect to to the entropy-entropy flux pair {rj,G) given by (1.8) if there 
exists a numerical entropy flux function Q{Wi, Wr) consistent with the exact entropy flux such that for 
any Wi , Wr we have 



G(Wr) + V^viWr) {J^r{Wl,Wr) - F{Wr)) < G{Wl, Wr) 

< G{Wi) + V^7j(Wi) {Ti{Wi,Wr) - F{wi)) . (A.l) 



In most situations, the modified hydrostatic reconstruction scheme presented here coincides with the 
original one. Thus, in most cases we will obtain an entropy satisfying scheme. 
To study the general case, we shall use a result presented in [7]: 
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Figure 24: Test 6: Number of cells needed to obtain -error < 0.008 
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(a) (b) 
Figure 25: Test 6: Number of cells needed to obtain -error < 0.008 
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(a) (b) 
Figure 26: Test 6; Number of cells needed to obtain L^-error < 0.008 
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Proposition A. 2. Consider a numerical scheme in the form 



with Ff^^j^ = T^iyVi, T^i+i) and 



j'+{Wi,Wr) = j'iwt,w;) + 





p{hr)-p{h;)+T+{Wl,Wr) 



(A.3) 



where J- stands for a given consistent numerical flux for the homogeneous shallow water system that 
verifies a semi-discrete entropy inequality for the entropy pair {r],G) in (1.8) and the interface values 
Wi,w* are derived from a local reconstruction procedure. Denote by F{w'[,w*) = {F^.F''). 
A sufficient condition for the scheme to he semi-discrete entropy satisfying for the entropy pair [rj, G) 
in (1.8) is that for some H* we have 



where 



£1 > 0, Sr< 0, (A.4) 



2 2 

+{ui-ut){F''-p{h^)) + uiT-, (A.5) 



£1 = jrh.i^g(^hi-h*i)-Hi+H* 
+{ui-ut){F''-p{h^)) + u, 
Sr = ■ (g(hr-h;-Hr + H*) + 



2 2 

+{ur - - piK)) + UrT+. (A.6) 

Consider now the modified hydrostatic reconstruction scheme proposed before and assume we are in 
Case 2. Assume now that F{wT'^J^^2^ ^4+1/2) ^ conservative flux for the homogeneous problem which 
satisfies a semi-discrete entropy inequality. Suppose for instance that hi — Hi + Hi+\ < 0. If we try to 
apply Proposition A. 2, we get that a sufficient condition for the modified hydrostatic reconstruction 
scheme to satisfy a semidiscrete entropy inequality is 

Q<F^- g{hi - h-^y^ -Hi + Hi+i) + «ii;+i/2 

= F'' ■ g{hi -Hi + Hi+i) - ghiUiihi -Hi + Hi+^). (A.7) 

If F defines a positive scheme for shallow water with flat bottom, then = J^''(0, ^^^^^2) ^ ^^'^ 
J^^ • gihi — Hi + -ffj+i) > 0. Unfortunately, the sign of —ghiUi{hi — Hi+ -ffj+i) will depend on the sign 

of Ui. 

As a consequence, we may say that the modified hydrostatic reconstruction is semi-discrete entropy 
satisfying in most cases, but eventually this property could be lost. 
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